{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>First full pipeline test on Antarctica</h1>\n",
    "<p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/proplot/config.py:1454: ProPlotWarning: Rebuilding font cache.\n"
     ]
    }
   ],
   "source": [
    "from icepyx import icesat2data as ipd\n",
    "import numpy as np\n",
    "import os\n",
    "import shutil\n",
    "import h5py\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import sys\n",
    "import pyproj\n",
    "import proplot as plot\n",
    "\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Earthdata Login password:  ··········\n"
     ]
    }
   ],
   "source": [
    "short_name = 'ATL06'\n",
    "spatial_extent = [31.5, -70.56, 33.73, -69.29]\n",
    "date_range = ['2020-03-30','2020-04-1']\n",
    "region_a = ipd.Icesat2Data(short_name, spatial_extent, date_range)\n",
    "\n",
    "username = \"JordiBN\"\n",
    "email = \"jordi.bolibar@univ-grenoble-alpes.fr\"\n",
    "\n",
    "region_a.earthdata_login(username,email)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#region_a.order_vars.avail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_a.order_vars.append(var_list=['count'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of data order requests is  1  for  1  granules.\n",
      "Data request  1  of  1  is submitting to NSIDC\n",
      "order ID:  5000000701699\n",
      "Initial status of your order request at NSIDC is:  processing\n",
      "Your order status is still  processing  at NSIDC. Please continue waiting... this may take a few moments.\n",
      "Your order is: complete\n",
      "Beginning download of zipped output...\n",
      "Data request 5000000701699 of  1  order(s) is downloaded.\n",
      "Download complete\n"
     ]
    }
   ],
   "source": [
    "region_a.download_granules('/home/jovyan/surface_classification/data') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "FILE_NAME = '/home/jovyan/data/processed_ATL06_20200330121520_00600712_003_01.h5'\n",
    "f = h5py.File(FILE_NAME, mode='r') \n",
    "\n",
    "count = f['gt1l/residual_histogram/count'][:] # has units of n_histograms, n_bins\n",
    "lat_mean = f['gt1l/residual_histogram/lat_mean'][:]\n",
    "lon_mean = f['gt1l/residual_histogram/lon_mean'][:]\n",
    "h_li = f['gt1l/land_ice_segments/h_li'][:]\n",
    "h_lat = f['gt1l/land_ice_segments/latitude'][:]\n",
    "h_lon = f['gt1l/land_ice_segments/longitude'][:]\n",
    "\n",
    "#latitude = f['/gt2r/heights/lat_ph']\n",
    "#longitude = f['/gt2r/heights/lon_ph']\n",
    "#height = f['gt2r/heights/h_ph']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Cropping the data far from surface in each histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = count[:, 200:550]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "775382cf22f340e0985523c5d679cd75",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:3: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    }
   ],
   "source": [
    "\n",
    "fig=plt.figure(figsize=(10,8))\n",
    "plt.title(\"Training data\")\n",
    "ax = fig.add_subplot(111)\n",
    "h = ax.imshow(np.transpose(data),vmin=0,vmax=30,cmap='inferno')\n",
    "plt.colorbar(h)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Plot research area of the above file</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** still needs track on this image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_root='/srv/tutorial-data/land_ice_applications/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cd ..; [ -d pointCollection ] || git clone https://www.github.com/smithB/pointCollection.git\n",
    "sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "import pointCollection as pc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Plotting track on map</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "07f0f7123d1b4d1eb3488a75dbe43340",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "spatial_extent_ps = [spatial_extent[0], spatial_extent[2], spatial_extent[1], spatial_extent[3]]\n",
    "\n",
    "## we will want to set colorbar parameters based on the chosen variable\n",
    "vmin=0\n",
    "vmax=6\n",
    "ticks=np.arange(vmin,vmax+1,1)\n",
    "\n",
    "plt.figure(figsize=(8,8), dpi= 90)\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo(central_longitude=0)) # choose polar sterographic for projection\n",
    "ax.coastlines(resolution='50m', color='black', linewidth=1)\n",
    "ax.set_extent(spatial_extent_ps, ccrs.PlateCarree())\n",
    "plt.plot(lon_mean,lat_mean,transform=ccrs.PlateCarree())\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot comparing mean_lon and mean_lon from histrograms with beam lat and lon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e222d87d785241ba89721fe9d74e792d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fd9e4400490>]"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(h_lon,h_lat,'ob' )\n",
    "plt.plot(lon_mean, lat_mean,'.r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Unsupervised learning of ATL06 residual histograms</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training data shape: (523, 350)\n",
      "\n",
      "Classified labels: [0 0 0 0 0 3 3 3 3 3 3 3 0 0 3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 3 0 3 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 3 3 3 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 3 3 3 0 0 0 3 3 3 3 3 3 3 3 3 0 0 3\n",
      " 2 3 3 3 3 2 2 2 2 3 3 3 3 3 3 2 2 2 2 2 2 3 3 3 3 2 2 2 1 1 2 2 2 1 1 1 1\n",
      " 2 3 3 2 2 1 1 1 1 1 1 1 2 3 1 1 1 3 3 3 3 3 1 1 2 2 3 3 2 2 2 2 1 1 1 1 2\n",
      " 2 2 1 2 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1\n",
      " 1 2 2 1 1 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 3 3 3 3 3 3 3 3\n",
      " 0 3 3 3 3 3 3 3 3 3 3 0 3 0 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 2 3 3 3 3 3 3 3\n",
      " 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 3 3 3 3 3 3 2 3 2 2 1 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1]\n",
      "\n",
      "K-means labels shape: (523,)\n"
     ]
    }
   ],
   "source": [
    "print(\"Training data shape: \" + str(data.shape))\n",
    "\n",
    "# Use int random_state in order to make centroid initialization deterministic\n",
    "kmeans = KMeans(n_clusters=4, random_state=1).fit(data)\n",
    "\n",
    "# Display classified labels\n",
    "print(\"\\nClassified labels: \" + str(kmeans.labels_))\n",
    "\n",
    "print(\"\\nK-means labels shape: \" + str(kmeans.labels_.shape))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the classified labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/proplot/ui.py:492: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  **kwargs\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2012130aa8c54231807a0c92fddfb621",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig1, ax1 = plot.subplots(ncols=1, nrows=2, share=0, width=5, height=6)\n",
    "\n",
    "fig1.suptitle(\"Classified labels along transect\")\n",
    "\n",
    "ax1[0].set_ylabel('Histogram frequency')\n",
    "\n",
    "ax1.format(\n",
    "        abc=True, abcloc='ul',\n",
    "        ygridminor=True,\n",
    "        ytickloc='both', yticklabelloc='left'\n",
    ")\n",
    "\n",
    "# Residual histograms\n",
    "ax1[0].imshow(np.transpose(data),vmin=0,vmax=30,cmap='inferno')\n",
    "ax1[0].colorbar(h)\n",
    "\n",
    "# Classified labels\n",
    "ax1[1].scatter(range(0,data.shape[0]), kmeans.labels_, c=kmeans.labels_, cmap='viridis')\n",
    "ax1[1].set_ylabel('Classification label')\n",
    "ax1[1].set_xlabel('Segments along track')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We display the labels on top of the raster map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-70.56 -69.29 -69.29 -70.56 -70.56]\n",
      "[33.73 33.73 31.5  31.5  33.73]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:14: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "43133013e344405eb337096ff52a22cf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([1114050., 1262050., 1773825., 1938825.]), 'origin': 'lower'}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Labeled transect')"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spatial_extent = np.array(spatial_extent)\n",
    "lat=spatial_extent[[1, 3, 3, 1, 1]]\n",
    "lon=spatial_extent[[2, 2, 0, 0, 2]]\n",
    "print(lat)\n",
    "print(lon)\n",
    "# project the coordinates to Antarctic polar stereographic\n",
    "xy=np.array(pyproj.Proj(3031)(lon, lat))\n",
    "# get the bounds of the projected coordinates \n",
    "XR=[np.nanmin(xy[0,:]), np.nanmax(xy[0,:])]\n",
    "YR=[np.nanmin(xy[1,:]), np.nanmax(xy[1,:])]\n",
    "MOA=pc.grid.data().from_geotif(os.path.join(data_root, 'MOA','moa_2009_1km.tif'), bounds=[XR, YR])\n",
    "\n",
    "# show the mosaic:\n",
    "plt.figure()\n",
    "MOA.show(cmap='gray', clim=[14000, 17000])\n",
    "ax.stock_img()\n",
    "plt.plot(xy[0,:], xy[1,:])\n",
    "# This still needs to be fixed in order to properly display the transect on the map\n",
    "x_polar, y_polar=np.array(pyproj.Proj(3031)(lon_mean, lat_mean))\n",
    "plt.scatter(x_polar, y_polar, c=kmeans.labels_)\n",
    "plt.title('Labeled transect')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
